{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from pathlib import Path\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "%run dataset.ipynb\n",
    "\n",
    "def select_gene_nonzeroratio(df, ratio):\n",
    "    nonzerocounts = np.count_nonzero(df.values, axis=0) / df.shape[0]\n",
    "    selected_genes = df.columns[nonzerocounts > ratio]\n",
    "    return selected_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasetname = 'onemillionv2'\n",
    "dataset = DATASET(datasetname)\n",
    "dataset.load_dataset()\n",
    "data_sc = dataset.data_sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "395\n"
     ]
    }
   ],
   "source": [
    "monocyte_ut = data_sc[(data_sc.obs['time']=='UT') & (data_sc.obs['cell_type_lowerres']=='monocyte')]\n",
    "monocyte_ut_df = pd.DataFrame(data=monocyte_ut.X.toarray(),\n",
    "                              index=monocyte_ut.obs.index,\n",
    "                              columns=monocyte_ut.var.index)\n",
    "mono_genes = select_gene_nonzeroratio(df=monocyte_ut_df, ratio=0.50)\n",
    "print(len(mono_genes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11482, 379) (194, 379)\n"
     ]
    }
   ],
   "source": [
    "bp_df = pd.read_csv('mono_gene_nor_combat_20151109.ProbesWithZeroVarianceRemoved.ProbesCentered.SamplesZTransformed.1PCAsOverSamplesRemoved.txt.gz',\n",
    "                   compression='gzip',\n",
    "                   sep='\\t', index_col=0)\n",
    "name_mapping_dic = pd.read_csv('features_v3_reformated_names.tsv',\n",
    "                          sep ='\\t',\n",
    "                          names=['geneid', 'genename']).set_index(['geneid'])['genename'].T.to_dict()\n",
    "\n",
    "bp_df['geneid'] = [item.split('.')[0] for item in bp_df.index]\n",
    "bp_df['genename'] = [name_mapping_dic.get(geneid) for geneid in bp_df['geneid']]\n",
    "bp_df = bp_df.dropna(subset=['genename'])\n",
    "bp_df = bp_df.drop('geneid', axis=1)\n",
    "bp_df = bp_df.set_index('genename')\n",
    "print(bp_df.shape)\n",
    "\n",
    "bp_trans_df = bp_df.T\n",
    "common_genes = list(set(mono_genes) & set(bp_trans_df.columns))\n",
    "selected_mono_df = monocyte_ut_df[common_genes]\n",
    "selected_bp_df = bp_trans_df[common_genes]\n",
    "print(selected_mono_df.shape, selected_bp_df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_mono_df.T.to_csv('sc_Expression.csv', sep=',')\n",
    "selected_bp_df.T.to_csv('bp_Expression.csv', sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create this fake pseudo time ordering because it's required to run the Beeline tool, but not used by GRNBoost2\n",
    "fake_timepoint_bp = pd.DataFrame(index=selected_bp_df.index)\n",
    "fake_timepoint_bp['time'] = np.arange(selected_bp_df.shape[0])\n",
    "fake_timepoint_bp.to_csv('bp_timepoint.fake.csv',\n",
    "                        sep=',')\n",
    "fake_timepoint_sc = pd.DataFrame(index=selected_mono_df.index)\n",
    "fake_timepoint_sc['time'] = np.arange(selected_mono_df.shape[0])\n",
    "fake_timepoint_sc.to_csv('sc_timepoint.fake.csv',\n",
    "                        sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform GRNBoost2 with BEELINE, see the yaml files in the same directory\n",
    "# python BLRunner.py --config config-files/config_bp_mono.yaml\n",
    "# python BLRunner.py --config config-files/config_sc_mono.yaml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Gene1_sc</th>\n",
       "      <th>Gene2_sc</th>\n",
       "      <th>EdgeWeight_sc</th>\n",
       "      <th>Gene1_bp</th>\n",
       "      <th>Gene2_bp</th>\n",
       "      <th>EdgeWeight_bp</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sorted_genepairs</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CCL3;CCL4</th>\n",
       "      <td>CCL3</td>\n",
       "      <td>CCL4</td>\n",
       "      <td>554.503642</td>\n",
       "      <td>CCL3</td>\n",
       "      <td>CCL4</td>\n",
       "      <td>55.157748</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CCL4;CCL3</th>\n",
       "      <td>CCL4</td>\n",
       "      <td>CCL3</td>\n",
       "      <td>480.484753</td>\n",
       "      <td>CCL4</td>\n",
       "      <td>CCL3</td>\n",
       "      <td>77.414467</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S100A9;S100A8</th>\n",
       "      <td>S100A9</td>\n",
       "      <td>S100A8</td>\n",
       "      <td>341.726427</td>\n",
       "      <td>S100A9</td>\n",
       "      <td>S100A8</td>\n",
       "      <td>104.542395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S100A8;S100A9</th>\n",
       "      <td>S100A8</td>\n",
       "      <td>S100A9</td>\n",
       "      <td>284.321568</td>\n",
       "      <td>S100A8</td>\n",
       "      <td>S100A9</td>\n",
       "      <td>65.915233</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S100A9;LYZ</th>\n",
       "      <td>S100A9</td>\n",
       "      <td>LYZ</td>\n",
       "      <td>221.872616</td>\n",
       "      <td>S100A9</td>\n",
       "      <td>LYZ</td>\n",
       "      <td>0.149064</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 Gene1_sc Gene2_sc  EdgeWeight_sc Gene1_bp Gene2_bp  \\\n",
       "sorted_genepairs                                                      \n",
       "CCL3;CCL4            CCL3     CCL4     554.503642     CCL3     CCL4   \n",
       "CCL4;CCL3            CCL4     CCL3     480.484753     CCL4     CCL3   \n",
       "S100A9;S100A8      S100A9   S100A8     341.726427   S100A9   S100A8   \n",
       "S100A8;S100A9      S100A8   S100A9     284.321568   S100A8   S100A9   \n",
       "S100A9;LYZ         S100A9      LYZ     221.872616   S100A9      LYZ   \n",
       "\n",
       "                  EdgeWeight_bp  \n",
       "sorted_genepairs                 \n",
       "CCL3;CCL4             55.157748  \n",
       "CCL4;CCL3             77.414467  \n",
       "S100A9;S100A8        104.542395  \n",
       "S100A8;S100A9         65.915233  \n",
       "S100A9;LYZ             0.149064  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sc_edges = pd.read_csv('sc_edges.csv', sep='\\t')\n",
    "sc_edges['sorted_genepairs'] = [';'.join(item) for item in sc_edges[['Gene1', 'Gene2']].values]\n",
    "bp_edges = pd.read_csv('bp_edges.csv', sep='\\t')\n",
    "bp_edges['sorted_genepairs'] = [';'.join(item) for item in bp_edges[['Gene1', 'Gene2']].values]\n",
    "\n",
    "sc_edges = sc_edges.set_index('sorted_genepairs')\n",
    "bp_edges = bp_edges.set_index('sorted_genepairs')\n",
    "concated_edges = pd.concat([sc_edges.add_suffix('_sc'), bp_edges.add_suffix('_bp')], axis=1)\n",
    "\n",
    "concated_edges.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SpearmanrResult(correlation=0.16937964029402044, pvalue=0.0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "concated_edges = concated_edges.dropna()\n",
    "spearmanr(concated_edges['EdgeWeight_sc'], concated_edges['EdgeWeight_bp'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Spearman r = 0.17')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAFNCAYAAACE8D3EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABDJElEQVR4nO3deXzU9Z348dc7d8h9AQlHwiUgiIgoWqyFStXSWm3toWu3uj3c7dpr3W5Xt/vbtdvt1nat23Pbdbuu7rZra6tttbUeKEhFqwIqECGEQMKREHKRgxyTTN6/P77fiZMwM5lJZiYz8H4+HmPm+/3OzPc9Ad987o+oKsYYY8KTMtUBGGNMMrGkaYwxEbCkaYwxEbCkaYwxEbCkaYwxEbCkaYwxEbCkaYwxEbCkaQAQkctE5EUR6RSRdhHZJiIXTXVcZwpxfENE2tzHN0VEgrw2Q0R+KSL1IqIism7M9d+LSI/fwyMiu+PxPYwlTQOISD7wW+B7QDEwC/gKMBDnONLieb/xRDmeW4HrgPOBFcB7gT8P8foXgI8Cx8deUNV3q2qu7wG8CPwiirGaECxpGoBzAFT1IVX1qmqfqj6tqrsAROQWt+T5Pbckuk9ErvC9WUQKROS/RKRJRI6JyD+LSKp7bYGIPOeWrlpF5KciUuj33noR+VsR2QWcEpGFbunqz0TkiIh0iMhfiMhFIrJLRE6KyPf93h/O53/RfW+niPxcRLIC/RL8vue/iUg7cFcUf8c3A99S1aOqegz4FnBLoBeqqkdVv62qLwDeUB8qIlXA24H/jWKsJgRLmgZgP+AVkQdF5N0iUhTgNWuAg0Ap8I/AoyJS7F57EBgCFgIXAFcCn3SvCfB1oAJYCszh9GR0I/AeoND9HN/9FgEfAb4NfBnYACwDPiwi74jg8z8MXA3Mwynl3RLid+H7ntOBr429KCJ/4ibuYI+5QT53GfCG3/Eb7rnJ+hjwB1U9FIXPMuFQVXvYA5yE8wBwFCdxPQbMcK/dAjQC4vf6V4A/BWbgVOOz/a7dCGwOcp/rgNf8juuBj/sdVwEKzPI71wZ8xO/4EeALEXz+R/2Ovwn8KMh7bwEOx+j36wWW+B0vcr+njPO+o8C6ENcPALdM9d+fs+mRUG1IZuqo6l7cEpiILAF+glPCu9F9yTF1/y91NeCU7iqBdKDJr18jBTjiftZ04Ls4Vcg891rHmNsfCRBSs9/zvgDHuRF8vn+7YK8bdzCBYomGHiDf7zgf6BnzO42IiFwGzAR+OcnYTASsem5Oo6r7cEqdy/1OzxrT2zsXp/R5BKekWaqqhe4jX1V9Vc+v45SoVqhqPk7nxthe48kstRXO50ciZCwictOYnuuxj2DV82qcTiCf891zk3Ez8Kiq9kzyc0wELGkaRGSJiPy1iMx2j+fglDD/6Pey6cDnRCRdRD6EU51/QlWbgKeBb4lIvoikuJ0zvjbHPJxS1kkRmQX8TZTDj/Xnj6KqP1W/nusAj8NB3vo/wO0iMktEKoC/xvmHKSARyfTrsMoQkSz/f7REJBv4UKjPMLFhSdMAdON0gLwsIqdwkuUenP+xfV7GaYdrxekg+aCqtrnXPgZkAG/iVI1/CZS7174CrAI6gd8Bj0Y59lh/frT8B/A4sBvnd/s79xwAIlItIjf5vb4GpxliFvCU+7zS7/p1ON95c0yjNqeRSTSpmLOEiNwCfFJVL5vqWIyZalbSNMaYCFjSNMaYCFj13BhjImAlTWOMiYAlTWOMiUBSzwgqLS3VqqqqqQ7DGHOG2bFjR6uqlgW6ltRJs6qqiu3bt091GMaYM4yINAS7ZtVzY4yJgCVNY4yJgCVNY4yJgCVNY4yJgCVNY4yJgCVNY4yJgCVNY4yJgCVNY4yJgCVNY4yJgCVNY2Kkz+Nle307fZ6QW5ebJGNJ05gYqW7s5Idb6qhu7JzqUEwUWdI0JkaWVRTw6XULWFZRMNWhmChK6gU7jElk2RmprK4qnuowTJRZSdMYYyJgSdMYYyJgSdMYYyJgSdMYYyIQs6QpIveLyAkR2eN3rlhEnhGRWvdnkd+1O0XkgIjUiMhVsYrLGGMmI5YlzQeAq8ecuwN4VlUXAc+6x4jIucANwDL3Pf8uIqkxjM0YYyYkZklTVbcC7WNOXws86D5/ELjO7/zPVHVAVQ8BB4CLYxWbMcZMVLzbNGeoahOA+3O6e34WcMTvdUfdc8YYk1ASpSNIApzTgC8UuVVEtovI9paWlhiHZYwxo8U7aTaLSDmA+/OEe/4oMMfvdbOBxkAfoKr3qepqVV1dVhZwW2JjjImZeCfNx4Cb3ec3A7/xO3+DiGSKyDxgEfBKnGMzxphxxWzuuYg8BKwDSkXkKPCPwN3AwyLyCeAw8CEAVa0WkYeBN4Eh4DZVtfW0jDEJJ2ZJU1VvDHLpiiCv/xrwtVjFY4wx0ZAoHUHGGJMULGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wELGkaY0wEgiZNEZkbz0CMMSYZhCpp/jpeQRhjTLIIlTQlblEYY0ySSAtxbZaIfDfYRVX9XAziMcaYhBYqafYBO2JxUxH5K+CTgAK7gT8DpgE/B6qAeuDDqtoRi/sbY8xEhUqabar6YLRvKCKzgM8B56pqn4g8DNwAnAs8q6p3i8gdwB3A30b7/sYYMxmh2jQ9MbxvGpAtImk4JcxG4FrAl6QfBK6L4f2NMWZCQpU0Pxxq2JGqHp7IDVX1mIjcAxzGaQJ4WlWfFpEZqtrkvqZJRKYHer+I3ArcCjB3ro2KMsbEV6ik+TucNkf/XnQFyoDpQOpEbigiRTilynnASeAXIvLRcN+vqvcB9wGsXr1aJxKDMcZMVNCkqarn+R+LSBVOG+MG4F8mcc8NwCFVbXE/91HgbUCziJS7pcxy4MQk7mGMMTEx7jRKEVkkIg8Av8fpTT9XVb83iXseBi4RkWkiIsAVwF7gMeBm9zU3A7+ZxD2MMSYmgpY0RWQ58GVgGfBN4BOq6p3sDVX1ZRH5JbATGAJew6lu5wIPi8gncBLrhyZ7L2OMiTZRDdwsKCJe4AhO2+ZpyTIRBrevXr1at2/fPtVhGGMC6PN4qW7sZFlFAdkZE+oCmTIiskNVVwe6Fqoj6BM4HT/GGBOx6sZOfriljk+vW8DqquKpDidqQnUEPRDHOIwxZ5hlFQV8et0CllUUTHUoURWqTfNxQpQ0VfV9MYnIGHNGyM5IPaNKmD6hquf3xC0KY4xJEqGq588HuyYia2MTjjHGJLZQ1fNU4MPALOBJVd0jIu8F/g7IBi6IT4jGGJM4QlXP/wuYA7wCfFdEGoBLgTtU9ddxiM0YYxJOqKS5GlihqsMikgW0AgtV9Xh8QjPGmMQTcmk4VR0GUNV+YL8lTGPM2S5USXOJiOxynwuwwD0WQFV1RcyjM8aYBBMqaS6NWxTGGJMkglbPVbUBp4f8Q8ASVW3wf8Qtwijq83jZXt9On2fS644kpDP9+xmTCIImTRH5d+CvgBLgqyLy/+IWVYz45sJWN3bG7B5Tmbji8f2MOduF6gi6HHinqt4JrOMM2LMnHnNhpzJxnalzfY1JJOP1nnsBVLWX0dteJCXfXNhYLlM1lYkrOyOVZRUFVDd2WhXdmBgJlTSXiMgu97Hb73i3X6+6GSOcxBzLKrxV0Y2JLes9nwKxXGfQqujGxFaopJmtqvsARCRTVQd8F0TkEiApe9ATQSwT25m6HJcxiSJU9fz//J6/NObav8cglrOGtT0ak7xCJU0J8jzQsYmQtT0ak5xCJU0N8jzQsYmQtT0ak5xCtWnOFpHv4pQqfc9xj2fFPLIznLU9GpOcQiXNv/F7PnafXNs31xhzVgq13cWD8QxkKiTzvszGmKkRyW6UirMQ8WZV/UmsA4uHM3Vf5lixf2SMiXw3ymLgoyKyXFXviFFMcZOsnTFTlbzsHxljQFQj6wh3N1zboaorYxJRBFavXq3bt59ZzavhJMTt9e1TkryspGnOFiKyQ1VXB7oWashRQL5FPJJRMqw3Gc74zakqIQebV58Mv1d/yRavSSyh1tMsDvBYICJfAarjGGPUJMOA8nASYjxWa4pEMvxe/SVbvCaxBK2ei8ghnM4f3+wfX0fQFuCrqtodjwBDibR6btXL2Ei232uyxWviL1T1PNSQo3khPvDnwEeiEFtc2YDy2Ei232uyxWsSS8Rtmq5LoxpFnFhbljFmsiaaNJOStWUZYyYr1OD2VcEuAemxCSe2xnayTLRty9rEjDl7hRrc/q0Q1/ZFO5CpMNHB2sk8yNsSvjGTE6ojaH08A4mHscku0vGOvoQzvzQ3KWcSQXInfGMSQcg2TRGpFJFS9/klIvJFEbkuLpHFwNgkGel4R1/COdjak1DjJCORrFNHjUkUocZp/gNwM874zJ8BG3DGaK4B3lDVL8QnxODiPY3SqrbGnB0mNE4TuAFnR8ppwGFgpqr2ikga8PokAyoEfgwsx0nKHwdqgJ8DVUA98GFV7ZjMfaLNxvcZY0JVz/tV1aOqJ4E6Ve0FUNUhwDPJ+34HeFJVlwDnA3uBO4BnVXUR8Kx7bIKwMafGTI1QSbNQRD4gItcD+e5z3/GEG8REJB+4HPgvAL/EfC3gW/j4QeC6id7jbGBjTo2ZGqGq588D17jPt/o99x1P1HygBfhvETkf2AF8Hpihqk0AqtokItMncY+AzqQ2yYl26JxJvwNjpkKoIUd/FsN7rgI+q6ovi8h3iKAqLiK3ArcCzJ07N6IbT+Vwm2gnq4m2r4b6HVhCNWZ8UzGN8ihwVFVfdo9/iZNEm0WkHMD9eSLQm1X1PlVdraqry8rKIrrxVA63SZTqdKjfQaLEaEwii3vSVNXjwBERWeyeugJ4E3gMZ4gT7s/fRPveU7kOZaKMjwz1O0iUGI1JZKHaNGPps8BPRSQDOAj8GU4Cf1hEPoEzxOlDUxRbTCTDcKVkiNGYqTZu0nT3BHoPzvjJkder6r0Tvamqvg4EGjh6xUQ/0xhj4iGckubjQD+wGxiObTjGGJPYwkmas1V1RcwjMUnHetvN2SicjqDfi8iVMY/EJB3rbTdno3CS5h+BX4lIn4h0iUi3iHTFOrB4OlOnJMb6e1lvuzkbhZM0v4WzJ9A0Vc1X1TxVzY9xXHE1XokpWZNqrEuCibaVsDHxEE7SrAX2aLA15M4A45WYkrUaaiVBY6Jf6AknaTYBW0TkThG53feIyt0TxHglpmgmn1B/gNH+w7WSoDHRL/SEkzQP4SzVlgHk+T3OCH0eL9tqW9lW2xI0WUUz+YT6A0zWEq0xiSzaNa6gK7ef9kKRPEBVtScqd46CaKzcvr2+na//fi8o3LlxacxnxIQaphPJEB4b7mNM7IRauX3ckqaILBeR14A9QLWI7BCRZdEOcqosqyjg9g2Luf1d58Sl7S9UqTWSEq2VSo2ZGuEMbr8PuF1VNwOIyDrgP4G3xS6s+MnOSGXtotKpDiNi1sljzNQIp00zx5cwAVR1C5ATs4hiKFmHDgVinTzGTI1wkuZBEfl/IlLlPv4ep3Mo6ViV1hgzWeEkzY8DZcCj7qMUZym3pFNekM380hzKC7KnOhRjTJIK2abpLgv3C1XdEKd4YmrzvmYefe0YlSXTuP7COVQ3djK/NJeDrT1nXS+09b4bMzEhS5qq6gV6ReSM6G1Yv2QGH7hgFuuXzBipqj+xu/GsrLJbU4UxExNO73k/sFtEngFO+U6q6udiFlWMNHX2cbD1FE2dfSO9z/NLc1lSnh+yF/pMLJVZ77sxExNO0vyd+0h680tzWb+4jPmluaO2dijODTyg3ZcsBwa93L+tfkp2sYwV29rCmIkJWj0XkWfdp+eq6oNjH3GKL6oOtvawuaaFg62nT2oKNBzJV4UFmXSpLBGHOyViTMYkulBtmuUi8g7gfSJygYis8n/EK8BoWlZRwE1rKtnR0MFze48HTJD+bXy+KuyqyqJJj4lMxDbERIzJmEQXdO65iHwQ+ARwGTB2greq6jtjHNu4JjL3/Ccv1fOvT9UwsyCLr73/vJEqaqzbLRNxXvmZ2FZrTDSEmns+7oIdIvL/VPWrMYlskiJNmn0eLy/VtVF7optF03O4dEFZQiaL7fXt/HBL3RnVhmpMMpnUgh2JmjAnorqxk5++3MCFlUVcuqCM6sZO+jzehGvbs55tYxJXODOCzhj+vef+7XnB2vYmk0wn816bV25M4jqrkqZ/7/myigI+vnYeXX2DdPUN8vG1VaeV7CbTUeJ7786GjoQqxRpjJiec9TT/N5xzycC/2pudkUpmego/2HKAH2w+QGZ66mklO19iHRj0Rpz0fPcCteRpzBkknMHtoxYcduejXxibcGJvYHCYnQ3trKosHlmAGDRg+6Evsf5wSx2Z6U6V2dfjPN6cdV8Vu8/jJTM9lYFBr3XuGHMGCJo0ReRO4O+AbL99zgXw4CxMnHR2NrTzT799k2npqXz5veeyuqp43AWIx3bK+Krd6xeXsbmmZdwkODZ5WueOMcktaNJU1a8DXxeRr6vqnXGMKWYGhpRhhSuWTh9V5Q41VnHsdMNI5qyH+hxjTHIat3quqneKyCyg0v/1qro1loHFQnf/ICd7PTS09rC94SSfSXeSZCTV5nDmrBtjzlzjJk0RuRu4AXgT8PViKJB0SfOlAy209HjYWtfGZ9cvYn5pLnubugL2nBtjTCDhdAS9H1isqgOxDibWZhRkk5EqXL6wdGQR4vu3HeLT6xbEbUykTV00JrmFtUcQkB7rQOJhSXk+M/OzuGzR9JEe8LEzb0INSo/GzKFEWyQj0WZDGZPoQi0N9z0R+S7QC7wuIv8hIt/1PeIXYvTkZ6VRlpfJ8a4+vv7EXvY2dZ4286a6sZPvbz7AIzuOnJZIQl0Ll39HUiySVaRJMNGSuDGJLlRJczuwA3gM+CrwonvseySdVZXF3LZ+EUfa+xjyDuOMoBptWUUBG5ZMZ9PeE6clklDXwuXrSDrY2hOTZBVpErR57sZEZtxVjhLZRFY5uvfpGn658wgfXDWH269cHLBdMVS7Y7TaJGPRttnn8bKzoQNQVlXGZu66tcmas8GkVjkSkd0ismvM4w8i8m8iUhL9cGNnZ0M7fzzYxiXzS7ll7TwAntt7gvv/cJD2Hs/I68ZbMMM3q2iqq9Zj+Tq2Ak0JjRarzpuzXTgdQb/H2SPoJvfxOPAH4DjwQMwiiwmhf8jLgeZu6lt7qG7s5Ou/38s9z+znu8/tD7hM3NjjnQ0d/NNvq7nn6f2TShyBks9kE1I8qtpWnTdnu3CS5lpVvVNVd7uPLwPvUNVvAFUTvbGIpIrIayLyW/e4WESeEZFa92fRRD87mKrSHBZOzyMjNYWa5m7ml+bywQtnUTQtndfqOwIuE3d6IlOmZaRxzYryCSUOXxIO1HM/2YQUjyXlbNk6c7YLJ2nmisga34GIXAzkuodDk7j354G9fsd3AM+q6iLgWfc4qjbva+blQ23MLpnGQ68c5uVDbSyvKOSr1y7nS+9ewrKKgtMS19jjVZXFfPk9S7lxTeVI4oikWu1Lwgdbe05LPmMTku9z23s8NizImAQRTtL8JPBjETkkIvXAj4FPiUgO8PWJ3FREZgPvcT/L51rAt8vlg8B1E/nsUC5dUMraBaUsKJ1Gc9cA2w60cM/TNdS3nhrpOBmbuMY7hsiGIvkn4fGSrS/BPrG7MWS13cZaGhM/4Wx38aqqngesBFaq6gpVfUVVT6nqwxO877eBLwHDfudmqGqTe88mYPoEPzuo3+9u5Kk3j9PU2ccXrljE2oWl9AwM8X8vH+ahVxrCatP06fN42VbbyrbaFuaX5o4MRRpv3Uz/pDteG6YvwW48ryJktd06Z4yJn1BLw31UVX8iIrePOQ+Aqt47kRuKyHuBE6q6Q0TWTeD9twK3AsydOzei9x472c/gkLJ1fyvXnD8HgDXzinjpYBu/ea2RjNQU5pXmcP+2+pEFPHwJyXfsG3IzMOjl3k37QeHOjUvZeF4FAANDQ9y/7UhYC4CM14YZzuIgfR4vA4PDNn/emDgJVdLMcX/mBXlM1FqcvdTrgZ8B7xSRnwDNIlIO4P48EejNqnqfqq5W1dVlZWUR3fiS+SUUZKeRl5lGfWsPP9hSS3VjNzetqeLalRVs2nsCkFGJzH9fIXirVAfC7RsWc/u7zmFZRcHIVhqZaWlx7V2OxzAjY8xbpnRwu1vS/KKqvldE/hVoU9W7ReQOoFhVvxTq/ZEObt9W28I3n6rh3PI8ZhVNY07RNNp6Brjugtkj1eWxg7bHbqcbbHD3RAZ9R2OrXhtsbkz0TXZw+zki8qyI7HGPV4jI30c7SOBu4F0iUgu8yz2OqqrSXGYWZPLUnia+/1wtT+5pYmttKwdbe8jOcFZV923r6zN2rrjv3NjXTWQoznjV83A6eGwIkDHxFU7v+X8CdwKDAKq6C2d9zUlT1S2q+l73eZuqXqGqi9yf7dG4h79fvXaUZ/ee4GSfl7QU4XBbLxuWTh9JWjsbOvj67/fyUl3bSLIKNFc8Wh0v4yU86+AxJvGEkzSnqeorY85NZnzmlNl19CRDw1CSk0pxbgYr5hRQUZjl9woFhfoAi2n4lwpjMSsmUKnSZt8Yk3jCWYS4VUScvWgBEfkg0BTTqGKkruUUACd7vWSmZ/B8bQvb6zu4632prF1UxqrKYu7cuJT5pbmsmFPI/NJcttW20tU/SNPJvpG2TyDq+/2M7aUH21fImEQUTtK8DWf3ySUicgw4hDMHPelcvrCUQy2neMc5JcwsnMa2ujaGgZrmblZVvpWcfMlqe307926qoblzgI5eD4jw8cuchT6i0QHj/xlWqjQmOYQzuP2gqm4AyoAlqnqZqjbEPrToWzm3iJK8DNLT03iupoX3rajgpovnsnV/60hbpf/MHt++6DdcPJuKwiyqSrJHPitYe+NEplRWN3aG3aFjs3+MmVrh9J7XichPgT8F5sQ+pNg50tHLyVMenqlupvFkP81dTpX78nPK6OrzUF6QzaKyXJ6qbh5JZKsqi7hgTjF3XbOMSxe8NS402HTISDpvJlK6jGbnkCVgYyIXTkfQucB/ACXAPSJyUER+FduwYmNwyIt3GIbcoal5WRk8sbuRx984xg821/FUdROv1LezsCxn1GD2QIPHg02HjCQRxmKYUiSsd96YyIWTNL04w428OHPFmwkyWyfRdfcNjUx2F6BnYJBHdh7jiiUzeM955ZQXZNHrGeKFA63sbXISSXlBNvNLcygvyA76uf6JLFbjJn2lQiBqn2/tqMZELpyk2YWzwMYh4GZVvVRV/zymUcXI9sMnAcjLSOFz71zAJfNL6PUMkZ4qbKtrIz8rgxsvriQ3Kw3f/kGb9zXzyGtHeWDboVHVWP+qbTwGmMeiWg7RS8DGnC3CSZo3AluBvwR+JiJfEZErYhtWbMwqynJ/ZnP+nELystLJzUpjbkkO6xeXsbQ8nxsvnssXr1wCKH0eLxvPq+D6C2az73h3VFdZD4d/YrZquTGJIZze89+o6t8Afw48AdwC/DbGccXEnMJpCLCv+RR//+tqak/0cPuGc8jPSmNzTcvIdEpQ7n1mPzsb2inOzeD2Kxfz+Q2LorrKejgm0rseDquWGzNx4fSePyIidcB3cFY++hgQ9a0o4qH2RDe+5UmG1cvvdjdxqPUUS8tH94TXNPcwrApIwPGYoRbtiGZvdKySm81XN2biwqme3w2co6pXqeo/q+rzqtof68Bi4bXDb01nL8+fxpVLnYWD/beeqG7sZGttC9evms3S8nwe2XGE7z93YGSBju317exsaA9YvfUvGUYjgVpyMybxhLty+xkxkK94WubI892Nnew61sXHLp07aqzl/NJcPrN+IRvPq+CJ3Y089WbzyKIe/mtpBioB+pcMk7nd0MZvGhNcOCXNM8b0QidpFmSmkJ0mbK1ppr6tD2CkROkrdR5s7WHTvhNcde4Mrr9wzsjScZ9et4BVlUUBS4D+JcNkbjdM5oRvTKydVUnzvIpCUoCUFPAiZGekUl6QSXVjJ5v2nRi1TJxv35+N51WE3FQtmGSuWidzwjcm1sZdsEOcTYFuAuar6j+JyFxgZoDl4hLeyd4BhoFUgV7PMDm5GeRnZbCsooDPrF84kiS217czMOhlc00LS8rzg+7Pc6ay1ZWMCS6ckua/A5fijNcE6AZ+ELOIYuiZvc0AtPYOc1FlEV+6ajHuincjW1nc+3QN39lUS7B2Sxjd5hfsuTHmzBRO0lyjqrcB/QCq2gFkxDSqGMlIfauq3NXvYWttKz96/uBI290Tuxt55LWjLJmZF7TdEka3+QV7PlHRTryWyI2JrnCS5qCIOCO+AREpY/R+5UljePitsLv7h3hyTxMz8zMYGPTS3uOhJDeTjcvLWT4rn50NbyWaYx193PVYNb/f3Xja7JxwV3QPdz/1aHfCWKeOMdEVTtL8LvArYLqIfA14AfiXmEYVI4U5bw05aukcYHAYGrv6uHfTfv7v5Xq+9rs32bSvmW8+VcM9T+8fSTQPbDvET16q51+e2Hfa7Jxgz8cam7yCJbNod8JM9POshGpMYON2BKnqT0VkB3AFzioW16nq3phHFgPFORnQ0gu4u8QBJ3sGOdk/RH1rD+2nPKSnCiU5mVyzonwk0dyydh6Dw8qaeUWjko9vZtD80lwOtvaMXAs0W2hs8gqWzAJ1wkxmlfiJduoE2n7DGBNe73kxzlJwD/mdS1fVweDvSkx7m7pGHednptB6aoDOPi+zinK4/sJZvHCgjY9cNIcb11QCsK22FVD+9uolo6ZR7mzooOZ4F1trW9mwdDqba1r4+Np5HGrtYdPeE3zmnQtHJZvJ9EhPRQKzYUfGBBbOHkE7cVZs78ApaRYCTSJyAviUqu6IXXjRNT0vne6BAcD5IlctKweG2Xmkk/PnFJCZlsabTd0srygkOyN1ZI8gFO7cuHQkYVU3dnLvphqGh+H6VbPYeF4F80pz2dPYyRO7mrj2gopxk024ibDP42VgcJib1sxlYNA7shRdrNmwI2MCCydpPgn8SlWfAhCRK4GrgYdxhiOtiV140dXc6Rl5rsCmNxupKMzlnYtnUFmSy/P7W7ht3QKWluePTKm8fYMzLGl+aS7b69tHOnxuW7eI+tYeNp5XQXFuBpnpKTy5p4n+IS+LZ+SNm9jml+ayfnHZyArxwfhWjl+/uIzNNS1kpk88mUVjMzhjznbhJM3VqvoXvgNVfVpE/kVVbxeRzFBvTDR9gzrquKNfye0fYkdDB282dfFmYyefveIcGk86c879p1Bur3cW6fj42nlkpqeQmeYsXLxiTiHFucUsqygYWYfTf2dLCJysDrb2hDV43ldNnl+ay5Ly/IBtquEmQWunNGbywuk9bxeRvxWRSvfxJaDDHYaUVEOPCnNOTyxHTvYx6B2murGTxTPzaTp5iu8/V8up/iGe2H18pHfbl7xA+eGWOgaGdFRJMTsjlbWLSlm7qOy0BBaop9w/GW6rbWVbbUvAnmpfNbk4N+O0nvlIhxNZO6UxkxdO0vwTYDbwa+A3wFz3XCrw4ZhFFgOpKaNLmgIsnpHDDRfN5f0rZ9HnGeJXrzdxvNvDG0dPcs7M3JEE40teqyqL3eQ5zK9eOzayl1Ao/gly7BYZB1t7uHdTDfc+sz/isZSRJsFkng9vTKIIZ2m4VlX9rKpeoKorVfUzqtqiqh5VPRCPIKOld+CtgnEKcOm8IoaGlV/uOMLbz5nOtRfM5v0rK8jPSqU0L4PLF5UCjBqv6Es8mWlpIDAwpOOOZ/RPkIFKnLdvWMzt7zpnJPmFO0bSkqAx8Rc0aYrI4yLyWLBHPIOMlp63+oEYBrr6h8hKTWX57AJqm7vYur+F8sJs5pXmsnFZOSvnFI9UgXc2dIxKZKsqi7jz3UvJTJOAVeRAiS9QyTBQtT7e+w8ZY8IXqqR5D/AtnF0o+4D/dB89wJ7YhxZ7e5q6mVsyjT3HOnmy2lls+Kpl5VxcVUzNiR5+/foxuvoG+fjaKnxtmb5E5lszE2TUcCCfQKu4AyMdMKESVrz3HzLGhC9o77mqPg8gIl9V1cv9Lj0uIltjHlkcpAIdpzwsKM1hzfxSuvsH2XW0g+0NHaycXcDjbxwjRYQ7Ny5lWUUBmempIx03vtWR7t9WP2o4kG/V9vmluaet4u7rtR6vF3vsGMlYDBWKRWK2IU3mbBDOkKMyEZmvqgcBRGQeUBbbsOLDC/yxvoO0wx10DXh5+WAb71tZTkoKrJxbhIqMTJ30JTL/Ae+3v2sxn163gPKCbMAZexkoIYY7hTKYWAwVisXgdRvSZM4G4STNvwK2iMhB97gKZzvfM0JGqjMnfenMPA61nmJGfjZXLytnT2MXj+48infYS35WOqsqnQ6X+aW5vOe8CqpKsllVWTQyhtM35jJYu6UviUykNBYsySZayc6GNJmzQTi9508Ci4DPu4/FvtlBySZ3zCqguenwycvm8Y/XLCc9LYWVs/N4ZOdR9jR2Mrc4m4LsdF6saxs1HOhgaw/bDrSSn50xkqj8k8V4PdoTaUsM9pmJ1i5pvfnmbBCq9/xLfofvU9U33MeAiCTl0nD+vefgLBVXMC2TVw62cd/zdWze38bRjj4efLGezLQ07v7ACr545RKuXDaTrj7PyFqaH187b1THTzjJwn+3S/891ifTg20lO2PiL1RJ8wa/53eOuXZ1DGKJu4GhYR577RjP17ZwyYJS3r+ygoKsNBQFhlm7qJT8rDQe3XmU7z17gJ0NHVQ3djIw5OXeTfvZ2dAx8lnjJUBfqXDsHuuTKSlayc6Y+AuVNCXI80DHSamlZ5Dp+Rl09w9ypO0UR0/2kZedTm5mujN4HQBhWkYa16ysYGBoiK/+9k1+9uphhgaH8fWgw/hV5UClwniUFG08pjHRFSppapDngY6TUmYq5Gan0+MZxovyh/0tXLawlLuuWcbS8ny21bbS1T/INSvKuW7lbDLT0ug45eGlA23MLc1haflbyW68VYsClQrjUVJMtHZPY5JdqN7z80WkC6dUme0+xz3OmugNRWQO8D/ATJyJOfep6nfcxY5/jtM7Xw982N3ELSZSAK8X6k/0kpOewoBnmKKcTNYvns7aRWVsq23ln35bjderpKYK5YXZNHX289dXLubV+nZ2Hz3J3qYu1rpTLcNdtcinvcfDE7sbR5aWC0c0e96NMRMTtKSpqqmqmq+qeaqa5j73HadP4p5DwF+r6lLgEuA2ETkXuAN4VlUXAc+6xzGT4gZS3dSFApcuKOEr7zuXSxc4SXBgaIhhhQvm5jOsyou1LXx/cy17jnWydmEJvYPD7Dl2kvYez2kdPDB+tfiJ3Y3826ZantjdGHbMgUqN493H2j2Nia5wxmlGlao2AU3u824R2QvMAq4F1rkvexDYAvxtrOIYcn9OyxBUlS37T5CeKlSW5PJSXSt5We6/C+L8u7J6XglZGWnsO95NeWE2KQJPVjeTnpbC73Y1cfu7zmHtorfG/Aca6O1fUtx4XgXAyM9wBCo12oByY+Ir7knTn4hUARcALwMz3ISKqjaJyPR4xLCgNI/mnn5auj38YsdRTvYN8WJdG9ecX06KwM7DHaRICm09A3x63UJeP3KSN5s6ueb8ChaU5XC4vY9hVcb2jYWT4D56adWk47fqtzHxNWVJU0RygUeAL6hql0h4HfIicitwK8DcuXMnHceRzj6uWFLGK4c6mJ6XxS1vq2L5rAIGvcN8YcMiMtNSqW89xeb9LWSkpeAZGuY/tx5kZkEWf3pJJS8dbOP6VbNZVVkUcHfKsTtS+o/xjLTKHKhU6Vs4JJFmBiXaTCVjoimcRYijTkTScRLmT1X1Ufd0s4iUu9fLcXbAPI2q3qeqq1V1dVnZ5KbA56RDT+8glSU5XHnuTA60dPPknuOkpwj/va2eth4P71w6g6uWl1OYnc7ju5oYHFb+cv1C/vbqxaxfMoMNS6az8byKkXGX3998gB9uOcD3nztwWo91dkYqmekp3L+tfkK92cFKlYnWQ55o8RgTTXFPmuIUKf8L2Kuq9/pdegy42X1+M84q8TF1atBpsuz1DHHOjDwWluXyYq2zpuZt6xdSUZhFe4+HB7Yd4oXaFtJShOf3t5CTkcqlC8po6uxjc00LB1t7ACepbVgynX3N3WxYOj1glTnQKu6TXXQ40aroiRaPMdEkqvEdcikilwF/AHbz1h5Df4fTrvkwznYah4EPqWp7qM9avXq1bt++Pex7V93xu4DnS3LSyUxPpa17gNzsdG59+zzmFufw7U37+cCqWWw70EZBdjot3QOU5WVysneQz29YFLBa7F81BYJW130btTlbZ2CdOcYkEBHZoaqrA12bit7zFwg+o+iKeMbiI+plXnEe3X0ehr1enth1nAvmFtLY2ceR9j7+/PL5LC0v4Ind7i6Vy2aMJMCxidNXGuzzeHlkxxE27T3BhqXT2VzTMpIUfXuZf3xt1UhytZKZMclhSnvPE0W3RznY1oNnSEnPTKE0L5NlFQUMDQ/zxuGTzC/NYVVlMddfOId5pbkEmj45tpRY3djJpn1Owtx4XsWo7Xd9e5l/et2CkZKnlTCNSQ5T0hGUaHLS4cOr5/Kh1bPJSEvhhdoWfrS1jrTUVE4NDnH/tkO8VNcCQE1zN/c8XTOyWIevR7yrb3DUNrzzS3NHOonGbr873pRLY0zisqQJdPYrh9tPsXZhKR+9pIqZBVmkp6Sw68hJSvMyaO3u52evHOGluhZ+8/oxWns8DAwNjbz/UGsP33uu9rR1N/07ifwFuxatxTVskQ5jYseSJuBVeOy1Jr69qZZzywuYUZDFvuZu+jxDdPYOIpLCC7UnqG/rY+N55WSlp4xs3buzoZ1N+05wzfmzRm3DG6oHeew1X5Lb2dAelaE6wXbQNMZMniVNV2l+BlcsKWNgaIij7X2kCrT3DjAtPY2h4WHetrCMq5bNJD1VyM1Io+lkHz/cUgcIn1m/kBsvnjtqG17/TqKxSWvs0CFfkgOJSoeQLymP3UHTGDN51hHkau7y8PCOYxxoOUVBVhqzCotoPdVPc3c/KSkpLJtVwH/+oY5XD7Zz8bxirlpezuKZ+QwMefHvGPIfchTuvPCxJc9Qw5jCmWHj34Pv2yHTGBMdVtJ0zS+ZxoBnkK01J5hTkkPPwBCluVl86rIF/MlFcwDl4VePcKyzl1cbOti8rxlQfrDlwKi2TP/ZMIGq6IHaG/1LnoFm00x0hk2gwfDW3mnM5FjSdPUMDNDtUdJThZxMofFkHzvqO3iloY1dxzpBISczjcLsDC6qLGLT3hOAcPuGxUHbMgMlrfHaGwO1dw4MDnPTmrmj9iWaKJviaMzkWPXcdaLHSUa9HuXFujYKp6Vxzow8ZhdlUXO8mz2N3czMy+L9F87iupWzAy7IAePvJ+5LigOD3rCq7r4xnesXl7G5poXM9NCfP15V3qY4GjM5ljTH8AKosHhmAXUt3Tz7ZgvtvYOc6Grl7edMp6pk2riJMZTx2hvHtoP6z1X3HyAfzHjtqJOJ3RhjSRNw5nSWF2TQ2Ons8Xve7ELqWrrxDCmIMi1TWDy9iNePdNDY3ktmWhqZ6SmUF2SzeV9zyC0rAs1FDzYFE04vCfonuXC20bCSpDGxZW2aOH3fvoQ5pyibOYWZ9Hq8qA5zpL2f3oFhZhXlcGrAywVVRfiG8jyw7dC4W1b4tyGObU8M1L44th000o4b297CmNiykqYrIwVQONHZx2O7jtPZN8iM/CwumVfM9atmkZGWSnqK8Mm3O6sSzS/N4boLZlFZMo31S2awvb49YDui/5TJ7IzUkVJgoEU7ArHtLIxJLFbSxKme52elkpmRyuqqIpZVFJCWIly+qITPvnMhFYXZ/M8f69l1rJP61lM8sO0Qj+w8ymuHO7j+wjls3tfMd56t5ZEdR04rEfpPmRw7tOj+bYfITE8NWSr0r277lzpt6JAxU8OSpqtwWibdA15ere/gRNcAN1xcyfT8LO56vJrWHg/vOa+C29Y7s2xqmru5ftVsLl1Qyr1P1/DEnuMsmZHHpr0nThvKE6yNMZy2x7E94aGq+saY+LCkidOm2dDWC4BnGIaGh1kzr4jf7TrOkfY+tje0se1AK/nZGayqLOZzVyzi9isX81JdK4+8dpRzZuRx6YISPvX2eaclwWADzMOZ4TM2MfonWuvwMWZqWNJ0DfotYF/feooHX2xg/eJSinMyyMtKZ+2CEsoLsqlu7KS8IJuHXjlMXlY6n1m3kDXzivjecweoOd4d1r3CLSUG60m3Th5jpo51BLlSgOJpqXT1eZlVnMUrh9rp6PUwMDjEIzuPkp+VDgLbDrQxvzSHh7cfYWZBFv/w3mXUNHfTMzDI47uaWDGncFSHTaBSZbilxFBjKq2DyJipYUnTNQzkZWXQ2tvHgrI8egeHqWnuIQVIT4Ws9BSqSnJYPCOf1p4BTvZ6SBPo6vewtbaFGy+uZPGM3HEHq0N0Bphb9dyYqWFJ009Dex8AzSf7ONHlGTk/4IVzywu4dEEp1Y2dPPDiIY6d7KPXM0RmZjqfumweqyoDV5vH2+t8onuE28weY6aGtWn68W2NOaRKuvubyclI5fKFJRRmptF4so+uviGWzMwlMy2FOUXT2NvYFXLY0Hh7nVsvuDHJxUqaAextPkVRdipVeVmkiuDxKg/vPMbhk/30eobo7h+iu99LxykPN15cNmpFImfvIB0peY43iH2y1eyJllSNMRNjJc0guvq9tPUOMLMgk8bOXlbMLuBIxykWTs9ldWUhswqyuO6C2ayqLOSlulae29vMQy83cM/T+05bX/P+bYcAGVnF3X9geqge8XAGsFtJ1Zj4spJmEF6Ftp4hXm3owDsMM/KyaO32sO1AKzmZ6Vy7soK27gF+tKWO3kEv/YPDFE/LYOOKctJTZGSnyUBLwQFhLws39nVjS5bWIWRMfFlJM4DsdBl5npuZxiXzSjjc3ssFcwpZOD2PoWEvzZ19/PqNRjIzUjl/VgGFWelce0EFyyvy2VbXNrLTpC+xgYxU0cNNdP7LwvlKnGNLljZ205j4sqQ5RjrQN6hMz0nliiVlLKsoJDcrhRPdHl6ub2dbbSsej7PR2mULS9l7rJM3jnXx/gtns3hGHkvLT98Hfew8c/+e72DV7/YeD4/sOML80lwOtvaE3ELDGBM/Vj0fY9D9eeKUl5fqWukbVMoLM0kTGHS718vysvAMebl8URl9g16WledRXpDJN5+qYdWcQt5+Tik/2FIHCnduXBpwC4vqxk4GBr3cv61+pPrtX/V+Yncj927aT0NbL59et/C0LTSMMVPDkmYIXq8iQOPJATJShfRhJTUFqps6+caTNQx6h+nqG6Sh7RQF0zLpOOXhlzuPMqckh9s3LAY0YInQV8W+aU3lyLJx/uc/vW4BG8+roKGtl33HuznY2mOJ0pgEYUkzhIFhZ9m4nHThlDs5PTNVKMrOYM28IrbWtlKYk84HL5zNdStnUVWSQ31rD9etnDWyknufx8tDrzTw+OuNfPaKReRnpzO/NJdPr1tAV5+HX712jHmlOaxdVHbapmy3X7l41KrvxpipZ22aYTg1qKTgJND01FTOm1XI3JJpDHqVdy8r52OXziM7I5X87DRuXFM5auuL6sZOfvN6I62nPNSe6OaHW+pGSo6ZaWnOhzr/Oa1TZ+xxe4+Hn7xUT3uPB2PM1LCkGUSK+/AVxQW4cG4hKQJ/2H+CE10eblu/kDXzi9nZ0M7OhvaA4yWXVRRw7coKSnMyWDQ9b9SUylWVRdz57qWsqiwKK6YndjeOu72GMSa2rHoehG9KpQDT0qAsP4u0VCEtVZiWkcor9W1UlU7j3mf2M6zw4dVzuGnNXFq6B/iP5+tYND2HlXOKOdjaw3UrZ7NiduHIRmo/3FI3shVvJG2VG8+rGPXTGBN/ljTH4QUGhqChvZ+G9n5KcjJ455IyXqhr46W6dgaGnM6gx99o5AOrZnH/C/W0dPczu3gab19UxhtHOvjilUtYu6gUmNy0yeLcDD56aVV0v6AxJiJWPQ+DF+cXlZMOl59TAiit3QP8sa6FvMw0Br1eVs4pZNCrZKQJJXmZvG1eMTvq22nt8TAwNDTq8wYGh9nZMPH9fcLdH8j2ETIm+ixphmkY8AzBK/XtvHiwg6KcDMrys+gfHOLyc2ZwUVURT+45zkcumss/X7uc9UtncNXymWSlpzAwpKNm9Ny7qYZvPlXDvU/XhNWpMzb5hTvf3OalGxN9Vj2PwKBCCoIIpKQJRzr6yUoThoehNDcDUJZXFJCZnsL3Nx+gLCeDjl4Pz+09zvFuD1edO4P1S2bwnvMqOHSim1/uPALA7VcuDmuvIN8g+EinYUZryJKtqGSMlTQjkp+dSr9nkI5eD8c6BkgB+oeUxpN9bNnXzJKZ+bx2pIOW7n4uX1jKi3WttPcMsq2ufWS3yqeqj/O73Y2snlfCZQvLeP3IyYAlTv/SZaCFjMOZbx7teelWcjXGkmZEuvq8tJzyMjQ0zDBv9bAPDENL9wCba47zvedq+epv3+Rg6yk+edkCLqws4rZ1C7h0QQkfu7SSQa+XYa/S1jNA98AQMwuyeGTn0dOGEfknKCfpKfc8vZ+HXm6YsjZKm/dujFXPJ6QvQM7q90JF/jS6+7rp6R3goVcOs7w8j0NtvdS25PPMvhOU52fxfG0Lf3LRHK5aXk5GWgoXzC2iJDeTvKw0ttW2sLS8gIOtPZQXZI+aYglCr2co4OZt8WLz3o1JwKQpIlcD3wFSgR+r6t1THFLY3jzeSf8QlOWk0zU4yP6WHoqmZVKam0l5QTY/3lpH66lBqpu6Ofd4J7967RieoWF2HzvJtgOtpAh85KI5bDvQxvrFZWyuaWFJeT7ZGQUMDA3zgVWzWTQ957SFP/zbGMNpd7S2SWMmLqGq5yKSCvwAeDdwLnCjiJw7tVGFr38IUgWqSqaRky5kpyp9nkGe3NXEH+taGfIOc+6MXHY2tPP73ccZ8iqDw8Pctm4hH7xwFtPSU6kqyeHT6xawfsmMkZLmzoZ2vvHkXp7Y1UTjyf6R+/lX4X1toMFmJvmztkljJi7RSpoXAwdU9SCAiPwMuBZ4c0qjioBX4dXDTjI6NQjgpaOvG453A3CyvwcFfr3zGJWl03h0x1FWzC4ElCuWTmdGfhYPbDtEbXM3fYNe6k70sHJuId5hpTQvg1+8eoSS3EzaegZYv2TGSBvjzoYO7nm6hquXzwi6H5HP2MWNfWt2TkXJ00q9JtkkWtKcBRzxOz4KrPF/gYjcCtwKMHfu3PhFFiXq/hwEDrT2kpeRysHWY3iHYVZRNjXNPfx2VxMKZKel0Nh5FBUhNUVoaO0lNVV4+WAbj+9qAvCbIaT0eoZ4cvdxLgyynbCPr21ye71TKvU1BYy3/UYsBNrSw5hElmhJUwKc01EHqvcB9wGsXr1aA7w+oZTkQP+AMzA+XSBnmgApLJtVRGlOBmvml3K0o5fmrn7esbiMc2bkk5WWQkuPh3ctncHAkJerlpePPM9MS6GqNJeF03NHzUFfVVnMP7x3GcHW8AzEv8S5pDx/SnrFrUfeJBtRTZy8IyKXAnep6lXu8Z0Aqvr1QK9fvXq1bt++PY4RGmPOBiKyQ1VXB7qWUB1BwKvAIhGZJyIZwA3AY1MckzHGjEio6rmqDonIZ4CncIYc3a+q1VMcljHGjEiopAmgqk8AT0x1HMYYE0iiVc+NMSahWdI0xpgIWNI0xpgIWNI0xpgIWNI0xpgIWNI0xpgIWNI0xpgIJNQ0ykiJSAvQEOHbSoHWGIQTa8kYdzLGDBZ3vCVi3JWqWhboQlInzYkQke3B5pQmsmSMOxljBos73pItbqueG2NMBCxpGmNMBM7GpHnfVAcwQckYdzLGDBZ3vCVV3Gddm6YxxkzG2VjSNMaYCTtrkqaIXC0iNSJyQETumOp4/InI/SJyQkT2+J0rFpFnRKTW/Vnkd+1O93vUiMhVUxM1iMgcEdksIntFpFpEPp/osYtIloi8IiJvuDF/JdFj9iciqSLymoj81j1O+LhFpF5EdovI6yKyPVniDkpVz/gHzoLGdcB8IAN4Azh3quPyi+9yYBWwx+/cN4E73Od3AN9wn5/rxp8JzHO/V+oUxV0OrHKf5wH73fgSNnacfahy3efpwMvAJYkc85j4bwf+D/htEv09qQdKx5xL+LiDPc6WkubI1sCq6gF8WwMnBFXdCrSPOX0t8KD7/EHgOr/zP1PVAVU9BBzA+X5xp6pNqrrTfd4N7MXZUTRhY1dHj3uY7j6UBI7ZR0RmA+8Bfux3OuHjDiJZ4z5rkmagrYFnTVEs4Zqhqk3gJCdguns+Ib+LiFQBF+CU3BI6dreK+zpwAnhGVRM+Zte3gS8Bw37nkiFuBZ4WkR3uFtyQHHEHlHDbXcTIuFsDJ5GE+y4ikgs8AnxBVbtEAoXovDTAubjHrqpeYKWIFAK/EpHlIV6eEDGLyHuBE6q6Q0TWhfOWAOem6u/JWlVtFJHpwDMisi/EaxMp7oDOlpLmUWCO3/FsoHGKYglXs4iUA7g/T7jnE+q7iEg6TsL8qao+6p5OithV9SSwBbiaxI95LfA+EanHaV56p4j8hMSPG1VtdH+eAH6FU91O+LiDOVuSZjJuDfwYcLP7/GbgN37nbxCRTBGZBywCXpmC+BCnSPlfwF5VvdfvUsLGLiJlbgkTEckGNgD7EjlmAFW9U1Vnq2oVzt/f51T1oyR43CKSIyJ5vufAlcAeEjzukKa6JypeD2AjTu9uHfDlqY5nTGwPAU3AIM6/tJ8ASoBngVr3Z7Hf67/sfo8a4N1TGPdlOFWnXcDr7mNjIscOrABec2PeA/yDez5hYw7wHdbxVu95QseNM2LlDfdR7ft/L9HjDvWwGUHGGBOBs6V6bowxUWFJ0xhjImBJ0xhjImBJ0xhjImBJ0xhjImBJ8wwkIl53RRnf47RVnURknW+lnKkiIn8hIh8b5zW3iMj3g1z7uxDv+5C7+tLmycYZLX6r/ewSkedFpNLvmorIt/yOvygid415/xsi8lAcQzYBWNI8M/Wp6kq/x91THVAgqvojVf2fSXxE0KSJM9b1L1V1vf9JEZnqqcPrVXUFzkykv/c7PwB8QERKA71JRJbi/P96uTtI3EwRS5pnEXHWFN0nIi8AH/A7X+auabhTRP5DRBp8//OKyEfFWX/ydfda6pjPvFhEHnWfXysifSKSIc66lQfd8wtE5El3wYY/iMgS9/xdIvJF9/lFbgnsJRH5V/FbWxSocN9fKyLfdF9/N5DtxvXTMTH9A87A+x+5n3WLiPxCRB7HWTiiWER+7d7vjyKywi+eB0XkabdU+AER+aZbOnzSnTI69nf6ORF50/2sn7nnckXkv/1KldcH+ON4idELUQzhbPvwV0H++P4E+F/gaeB949w/R5w1Wl8VZ+3Na93z2SLyM/e1PxeRl0UkaXaBTBhTPbreHtF/AF7emqHzOvARIAtn9ZhFOIsiPMxbs0q+D9zpPr8aZ5ZPKbAUeBxId6/9O/CxMfdKAw65z+/BmbK6FngH8JB7/llgkft8Dc4UQIC7gC+6z/cAb3Of3427tihwC3AQKHC/QwMwx73WE+J3sAVY7fcZR3FnnQDfA/7Rff5O4HW/eF7AWS7ufKAXd0YKzpzp6wLcpxHIdJ8Xuj+/AXzb7zVF7s963HUlcVYsutXvNT1AvvuaAuCLwF1+1/cDlTjTEB8b5/7/AnzUd859bw7OWpz3u+dX4CTq1VP99zXZHlNdVTGx0aeqK/1PiMhKnORW6x7/BPAt03UZ8H4AVX1SRDrc81cAFwKvirNyUTZvLayA+/ohcVbZXoqzEMO9OIsqpwJ/EGcFpLcBv5C3Vj/KHBNbIZCnqi+6p/4PeK/fS55V1U73tW/iJA//5cPC8Yyq+tYsvQy43o3/OREpEZEC99rvVXVQRHa73+FJ9/xuoCrA5+4CfioivwZ+7Z7bgDM/HPceHX6v3ywiM3B+j/7Vc9RZIep/gM8Bfb7zInIR0KKqDSJyFLhfRIrczw10/ytxFvf4onucBczF+XP5rnuvXSKyK8jvyoRg1fOzS7A5s8HWchPgQX2rbXSxqt4V4HV/AN6NM3d+E05SugzYivN37KSObmNdGub9fQb8nnuZ2JKGp8a5n+93MwCgqsPAoLrFMpw1LAPd9z3AD3D+cdnhtpkKwX/X63GSfjXwTwGufxunPda/3fJGYIk4KxzV4ZRIfVX+YPe/3u/3PVdV9475nmaCLGmePfYB80RkgXt8o9+1F4APA4jIlYBvv5ZngQ+Ksw6ib1+XSk63FfgC8JKqtuAsxrAEqFbVLuCQiHzI/QwRkfP93+yWmLpF5BL31A2EZzBQO2MYtgI3ufGsA1rdOCMiIik4TQWbcRYHLgRycdodP+P3uiL/96lqH87v62MiUjzmWjtO08kn/O7xIWCFqlaps8rRtcCNIe7/FPBZcYv2InJBgO+9HKeKbiJkSfPM5Osg8T3uVtV+nOr478TpCGrwe/1XgCtFZCdOibEJ6FbVN3GqkE+7VblncPYFGutlYAbO/5TgVBl3+ZXSbgI+ISK+lW4CbTXyCeA+EXkJp6TUGcb3vA/YNbYjKAx3Aavd73Q3by1RFqlU4CduVf414N/UWaPzn4EiEdnjfuf1Y9+ozmrlDwG3Bfjcb+G0KYNTpT6mqsf8rm/F2UtnVpD7fxWnXXaX26H2Vfd9PwRy3e/9JRJtybUkYascGUQkE/C67ZOXAj8c2yYahxhy1d27R5xxpeWq+vl4xnC2EZEtOB1x26c6lmRiHUEGnE6Ch93qngf41BTE8B4RuRPn72QDTo+3MQnHSprGGBMBa9M0xpgIWNI0xpgIWNI0xpgIWNI0xpgIWNI0xpgIWNI0xpgI/H/zV4rDA7ANLQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "coef, p = spearmanr(concated_edges['EdgeWeight_sc'], concated_edges['EdgeWeight_bp'])\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.scatter(concated_edges['EdgeWeight_sc'], concated_edges['EdgeWeight_bp'], s=1, alpha=0.5)\n",
    "plt.xlabel('Edge weight from scRNAseq')\n",
    "plt.ylabel('Edge weight from BLUEPRINT')\n",
    "plt.title(f'Spearman r = {coef:.2f}')\n",
    "# plt.savefig('grnboost2_sc_bp_comparison.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
